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Quasi-chemical theory is utilized to analyze the relative roles of solute polarization and size in 
determining the structure and thermodynamics of bulk anion hydration for the Hofmeister series 
CD, Br~, and D. Excellent agreement with experiment is obtained for whole salt hydration free 
energies using the polarizable AMOEBA force field. The total hydration free energies display a 
stronger dependence on ion size than on polarizability. The quasi-chemical approach exactly par- 
titions the solvation free energy into inner-shell, outer-shell packing, and outer-shell long-ranged 
contributions by means of a hard-sphere condition. The inner-shell contribution becomes slightly 
more favorable with increasing ion polarizability, indicating electrostriction of the nearby waters. 
Small conditioning radii, even well inside the first maximum of the ion-water(oxygen) radial distri- 
bution function, result in Gaussian behavior for the long-ranged contribution that dominates the 
ion hydration free energy. This in turn allows for a mean-field treatment of the long-ranged con- 
tribution, leading to a natural division into first-order electrostatic, induction, and van der Waals 
terms. The induction piece exhibits the strongest ion polarizability dependence, while the larger- 
magnitude first-order electrostatic piece yields an opposing but weaker polarizability dependence. 
The van der Waals piece is small and positive, and it displays a small ion specificity. The sum of 
the inner-shell, packing, and long-ranged van der Waals contributions exhibits little variation along 
the anion series for the chosen conditioning radii, targeting electrostatic effects (influenced by ion 
size) as the largest determinant of specificity. In addition, a structural analysis is performed to 
examine the solvation anisotropy around the anions. As opposed to the hydration free energies, the 
solvation anisotropy depends more on ion polarizability than on ion size: increased polarizability 
leads to increased anisotropy. The water dipole moments near the ion are similar in magnitude to 
bulk water, while the ion dipole moments are found to be significantly larger than those observed 
in quantum mechanical studies. Possible impacts of the observed over-polarization of the ions on 
simulated anion surface segregation are discussed. 

PACS numbers: 82.60.Lf,87.16.A-,61.20.Ja,64.70.qd,64.75.Bc 
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Specific ion or Hofmeister effects are exhibited in a wide range of physical, chemical, and biological systems^ often 
with significant impacts on bulk properties!^ Various explanations have been put forward as to the origins of these 
fascinating phenomena. Proposed contributing factors include: ion size^^ polarizabilityfSi^ dispersion forcesji and 
ion- water interactions^ While there are known specific ion effects observed for cation sequences^ a great deal of recent 
attention has focused on anions due to their large size, polarizability, and dispersion interactions^ As an example, 
Petrache et al£ presented experimental results that show a large increase in the swelling of biological membrane 
multilayers when the CD ion is replaced with Br~. The increase appears to be due to selective adsorption of Br~ 
ions relative to CD ions at the membrane surface. Other systems and phenomena impacted by specific ion effects^ 
include ion solvation free energies, ion activities, surface tension increments, bubble interactions, colloid interactions, 
pH measurements, buffers, zeta potentials, ion channels^ and protein stability^^iii 

The segregation of certain ions to the liquid-vapor interface of water provides another example of ion specificity^ 
Experiments suggest that larger, often more polarizable, anions exhibit enhanced occupancy at the interface relative 
to bulk concentrations i 13 i 14 i 15 ' 16 i 17 ' 18 These findings have implications for the understanding of atmospheric chemical 
reactions, sea-salt aerosols, and ions in biological systems. Molecular dynamics simulations have implicated the anion 
polarizability as a key determinant of surface specificity— d 3 i 19 i 20 i 21 The inclusion of anion polarizability for the halides 
appears to lead to an ion free energy minimum at water droplet or liquid-vapor interfaces in molecular dynamics 
simulations il&S^ Stuart and Berne^ 3 - utlizing a fluctuating-charge water model and a Drude oscillator polarizable 
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ion model, came to a different conclusion: the larger dipoles in water due to induction seem more important for ion 
surface affinity than polarization effects on the chloride ion itself. As discussed above, other theoretical work has 
targeted dispersion interactions,—^ ion-water interactions , 25 ! 26 and ion size.— £ There is clearly not a single agreed- 
upon mechanism that explains the observed experimental results. 

The above discussion begs the question of what level of theory is required to capture the essential physical effects 
leading to ion specificity. If ion size and/or polarizability are the major determinants, then classical models with 
accurate representations of polarizabilitics should be sufficient. The classical models assume, however, that fixed point 
multipolc expansions and induced point dipoles accurately capture the quantum mechanical charge distributions in 
the condensed phase. They also assume a pairwise form for the van der Waals dispersion interactions, thus neglecting 
possible many-body contributions. If complex charge rearrangements occur and/or there are appreciable many- 
body contributions to the dispersion energy, more accurate quantum mechanical treatments may be required. Ab 
initio molecular dynamics (AIMD) simulations at the gradient-corrected density functional theory (DFT) level are a 
common current choice for quantum modeling in the condensed phase. A well-known deficiency of DFT, however, is 
an improper handling of dispersion interactions. If those interactions play a signifcant role in ion specificity, a higher 
level of theory is required. As examples of the possible impact of a lack of proper dispersion interactions in DFT, 
the density of liquid water in AIMD modeling deviates from experiment by almost 20% (21 and cluster energetics are 
incorrectly estimated.— Addition of dispersion interactions at an approximate level improves the agreement of the 
DFT models with experiment^ 

An alternative approach is to measure or compute accurate frequency-dependent polarizabilitics for the ions and 
bulk water, and then to utilize those polarizabilitics in a Lifshitz approach for dispersion contributions to solvation 
free energies (or surface interaction free energies) 1 30 i 31 i 32 From these quantum chemical studies, it is apparent that 
a wide range of imaginary frequencies is required to obtain converged interaction energies. This strategy has sound 
fundamental underpinnings (as a many-body theory), but can fail to capture specific interactions of the ions with 
nearby waters in the condensed phase environment. Another avenue developed recently is to compute free energies 
with a classical model and then correct the classical free energies via a perturbation formula ; 33 ' 34 ' 35 the quantum 
correction is estimated from averages of the Boltzmann factor for interaction energy differences between quantum and 
classical models calculated on configurations from the classical simulation. Since quantum mechanics is not required to 
generate the trajectories, it is feasible to utilize correlated quantum mechanical methods (such as MP2-level theory). 
Mixed quantum mechanical/molecular mechanics (QM/MM) methods provide another viable approach for directly 
modeling ion solvation ! 36 ' 37 AIMD-DFT simulations coupled with thermodynamic integration methods have recently 
been employed to compute ion solvation free energies in water 

In this paper, we take a step back from the difficult task of a high-level quantum treatment. As a first step, we 
examine a well-tested classical polarizable model for ions in bulk water. The purposes of this study are: 1) to illustrate 
the utility of quasi-chemical theory in disentangling physical components involved in ion hydration, 2) to determine 
the major factors driving ion specificity for hydration free energies in a classical polarizable model, and to test that 
classical model's accuracy, 3) to examine local solvation structure as functions of ion size and polarizability in the 
classical model, and 4) to develop methods and perform benchmark classical calculations which set the stage for 
future, more accurate QM/MM-based calculations. By varying ion size and polarizability independently, we can gain 
insights into their relative contributions in ion hydration. This analysis is of course somewhat artificial physically, 
since the ion size and polarizability are typically related to each other. 

The model we choose for the present ion hydration studies is the AMOEBA force field. 39 ' 40 ' 41 This model includes 
fixed partial charges, dipoles, and quadrupoles on each water atom, as well as induced dipoles on each atom. The 
ions possess charges and induced dipoles. The van der Waals (vdW) interactions are handled with a buffered 14-7 
potential^ which has a replusivc wall intermediate in steepness between the Lennard- Jones and Buckingham forms; 
this potential yields a relatively accurate estimate of dispersion interactions at the pairwise level. The AMOEBA 
model produces a density maximum for pure water at 290 K and 1 atm. In addition, the model accurately reproduces 
the experimental density, heat of vaporization, radial distribution functions, magnetic shielding, self-diffusion, and 
static dielectric constants of bulk water. Also, ion parameters used with the AMOEBA water model have been shown 
to reproduce experimental hydration free energies^ 9 - Alternative polarizable force fields exist, including the Dang- 
Chang model^ POL3^- POL5/TZ^i NEMO, 7 . TTM2^ DPP^ Drude oscillator models^ and fluctuating charge 
models such as TIP4P/FQi 48 i 49 ' 50 Defusco, Schofield, and Jordan^ 6 - have recently provided a careful comparison of 
several polarizable water models by analyzing cluster structure and energetics. 

The AMOEBA force field includes a Thole-type damping function that avoids a 'polarization catastrophe' at short 
range. Masia, Probst, and Rey 5 -^ showed that inclusion of a well-fitted damping function is crucial to model ion-water 
dipolc moments properly By comparison with AIMD-DFT simulations, Masia 5 - 2 - further showed that simply reducing 
the ion polarizability in the condensed phase is not sufficient to represent the distribution of instantaneous dipole 
magnitudes accurately. In that work, a longer-ranged Gaussian form for the damping function was employed relative 
to the Thole version used in the AMOEBA force field, reflecting the diffuse electron densities of these ions^ 3 - The 
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FIG. 1: Schematic representation of the quasi-chemical approach to solvation. 



Masia model accurately reproduces the distribution of ion dipole magnitudes from AIMD simulations. Finally, in a 
more extensive AIMD study, Guardia, Skarmoutsos, and Masia£^ observed slightly decreased water dipole moments 
near halide ions relative to bulk, and they compared the anion dipole distributions to a recent classical polarizable 
simulation (without damping) by Wick and Xantheasi 5 - 5 - The results show an over-polarization of the ions in the 
undamped simulation. Other classical polarizable models have employed reduced ion polarization modelsi 5 - 6 - We will 
show below that the AMOEBA model also over-estimates the ion dipole magnitudes, even with the Thole damping 
function. This is likely due to the shorter range of the damping function relative to the model developed by Masia^ 
To address the aims stated above, here we utilize quasi-chemical theory to examine the influences of ion charge, 
size, and polarizability on ion hydration. We explore the relative contributions of first-order electrostatic, induction, 
and van der Waals interactions to the free energies, and also perform a structural analysis of the local solvation shell. 
The remainder of the paper is structured as follows. In the next section, wc briefly review quasi-chemical theory as 
applied to our analysis of ion hydration. Then we discuss the computational methods employed in the present study. 
The results of the simulations are presented, followed by a statement of our conclusions, discussion of the connections 
to related work, and future directions for investigation. 



II. QUASI-CHEMICAL THEORY 

The quasi-chemical theory (QCT) breaks down the solvation free energy into a three-step process . 57 ' 58 ! 59 ' 60 ' 61 ' 62 ! 63 
The steps include insertion of a hard sphere whose size is gauged by the length parameter A, addition of the solute in 
the created cavity, and finally removal of the hard sphere to allow close contact of the solute and solvent (Figure [1]). 
Each term depends on A, but the final result for the free energy is independent of A. 

In the first step, a molecular-sized cavity is created in the solvent to accommodate a hard sphere solute. The 
(outer-shell/hard-sphere, or packing) free energy required is given by 

^°o x s,HsW = -lnpo(A), (1) 

where po(X) is the probability of spontaneous cavity formation (i.e. the solvent center closest to a defined point in 
solution has a distance further than A). 

Next, the solute is introduced into the cavity center. This step can be considered in either the forward 
(uncoupled, 0— >coupled,l) or reverse (1— >0) directions, giving two equivalent formulas for the outer-shell, long-ranged 
free energy contribution: 

^os,lr(A) = - In (c-? AU \R min > A) (2) 
= ln(e^|i? min >A) r 

For A values beyond a critical value A m f , it is observed that the distribution of the interaction energies assumes a 
Gaussian form. 60 ' 63 In this regime the contributions to AU come from many interacting solvent molecules with similar 
interaction magnitudes. This situation makes the use of a quadratic perturbation formula a good approximation to 
the true Mos,lr : 

M) X s,lr(A) « «/3AC/|i? min > A) + (PAU\R min > A)J /2. (3) 
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This equation can be obtained by performing a cumulant expansion of the two forms of Eq.[2]to second order and then 
averaging. Since the distributions are assumed to be Gaussian, the fluctuation term is independent of the coupling 
state of the solute, and those terms cancel in the averaging. 

An important point to notice about the Gaussian approximation to ^Iqq lr is that it implies a linear change in the 
average solute-solvent interaction energy with respect to the solute coupling (traditionally A, but here we'll use 7): 

(AC/|A) 7 = (AC/|A) - / 3 7( r 2 , (4) 

with 

a 2 = (5A[/ 2 |A) 7 = ((AU\X) (AU\X) 7 ) //? 7 . (5) 
The mean-field estimate of the long-ranged contribution is then 

^os,lr(7; A) = ^ ((A[/|A> + (A£7|A> 7 ) . (6) 

An example of a Gaussian approach is the Born model of solvation, where the coupling constant 7 scales the ion 
charge, and the Born radius is determined by the fluctuation contribution. Below we will compare the electrostatic 
part of the long-ranged free energy contribution to the Born model prediction. 

In the final step, the hard sphere condition is eliminated, allowing the closest solvent molecules to form an inner 
solvation shell around the solute, yielding a chemical binding contribution to the free energy: 

(3^ = \nx (X), (7) 

where xq(X) is the probability of spontaneous cavity formation around the solute center with the solute included in 
the sampling. In practice, we use the Bayesian methods developed in Ref^ to compute po(X) and xq(X). In this 
approach, the observed solvent occupancies in shells around the cavity or solute center are employed to compute the 
free energies. An alternative direct approach for computing xq and po has been developed by Sabo et al.M- 

In a recent study, Merchant and Asthagiri^i provided an extensive classical analysis of ion solvation by carefully 
examining the thermodynamic contributions of particular coordination states to the inner-shell chemical term (Eq.[7]). 
This analysis showed that the inner-shell term exhibits significant stepwise increments up to the most probable 
coordination number, with much less variation beyond that number. The authors concluded that the ion influence 
on the solvent matrix is local based on these results. This important finding is consistent with the recent general 
discussion given by Collins, Neilson, and Enderby^ 

In the present paper, we adopt a strategy that to some extent is counter to the original goals of quasi-chemical 
theory 57 i 59 For the free energies computed here, we increase the value of A from until a stable total solvation free 
energy is obtained, and term this value A m f (where 'mf refers to mean- field). Further increases of A do not alter 
the total free energy, indicating that by this value of A the distribution of interaction energies is nearly Gaussian. 
We find that, for the spherical ion hydration problems examined here, a A m f value well inside the first maximum of 
the ion- water (oxygen) g(r) is sufficient to yield Gaussian behavior for the long-ranged contribution. Previous work 
has shown that non-Gaussian behavior is linked to high energy impacts of the waters with the solute, and a small 
level of conditioning removes the resulting high-energy tails i 60 i 63 For the systems studied here, a good choice for A m f 
is the average distance of the closest water oxygen to the ion. While we do examine the A dependence for a range 
of quantities, the above strategy was employed for most of the free energy calculations. We also note that such a 
small value of A m f is not accurate for a molecular solvation problem such as the water molecule in water^ For these 
anisotropic, molecular, hydrogen-bonded cases, the surrounding solvent must be pushed further out from the solute 
before Gaussian behavior is observed. 

With this choice for A m f, the contribution of /ijg is very small (roughly -0.5 kcal/mol) for all three anions (Cl~, 
Br~, and I - ) and exhibits little ion specificity. The packing contributions /^Qg HS are 5.3, 6.3, and 7.7 kcal/mol for 
the Cl~, Br~, and I ions, respectively (below). Those packing contributions depend only on the water model and 
A. So the /ijg + /ig'g HS combination for the chosen A m f values is positive and exhibits a small ion specificity mainly 
through the size-dependent packing term. The rest of the hydration free energy is then included in the outer-shell 
long-ranged contribution. With this choice for the spatial division, the outer-shell long-ranged piece includes, besides 
interactions with distant waters, contributions from waters close to the ions. 

A physically based interpretation of quasi-chemical theory often includes the first solvation shell as the inner 
shell, handles that spatial domain with accurate quantum calculations, and treats the longer-ranged contributions 
cither with a dielectric model or with classical simulations^ Here our aim is to start to disentangle the various 
energetic contributions to the ion hydration. Pushing most of the problem into the long-ranged contribution that 
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can be handled with a mean-field treatment allows for a division of the free energy into various parts: first-order 
electrostatic, induction, and vdW (defined below). There is no unique division of the free energy, but the present 
approach would appear to provide a natural and physical division. A side product is that, once the inner-shell and 
packing contributions have been determined, the mean-field long-ranged contribution can be computed in relatively 
short simulation times since it simply involves energetic averages for the coupled and uncoupled cases. This division 
should also prove useful in ongoing quantum mechanical studies of ion hydration. 

III. COMPUTATIONAL METHODS 

Molecular dynamics (MD) simulations were performed for several solutes in water using the NPT ensemble at 1 bar 
and 298.15 K. The Amber 10 package^ implementation of AMOEBA was employed for the simulations. All simula- 
tions, involving 215 waters and the various solutes, employed velocity Verlet integration with a timestep of 1.25 fs, an 
induced dipole convergence parameter of 1CP 4 Debye, a tapered cutoff ending at 9 A for the vdW interactions, and 
their associated long-range corrections. Each system was equilibrated for at least 60 ps using standard weak-coupling 
algorithms with temperature and pressure relaxation constants of 0.2 and 0.5 ps, respectively. This equilibration 
period is consistent with previous AMOEBA simulations of water and ions i 39 ' 40 i 66 We observed no significant drift in 
any of the free energy components considered in the next section during block averaging of the data. The correlation 
times for the distance of closest approach for the closest water to the ion are less than 0.5 ps for all systems, and the 
average lifetimes for water occupancy in the first solvation shell are on the order of 10-20 ps. In addition, we note 
that previous production runs were typically used as input for the next simulation. Since our simulations totalled 
over 20 ns in production runs, we are confident that the systems are properly equilibrated. Production runs of 1 ns 
duration used Andersen velocity randomization every 0.25 ps and a pressure relaxation constant of 1.5 ps during the 
collection of 2000 samples for analysis. Some additional simulations were performed using the nonpolarizablc SPC/E 
water model for comparison. 

The anionic solutes studied were classified first by their vdW parameters as: 1, CU-like, 2, Br~-like, or 3, I -like, 
and second by their polarizabilities (A 3 ): 0, a = 0.0 (nonpolarizable) , 1, a = 4.0 (Cl~-like), 2, a = 5.65 (Br~-like) or 
3, a = 7.25 (I - -like) . All ion simulations employed the vdW parameters for the individual ions from the AMOEBA 
code with no changes of well-depth or size. Of course, the dispersion component of the vdW interactions depends 
on the ion polarizability; here, we omit this dependence and defer a careful analysis of the dispersion interactions to 
ongoing quantum mechanical studies. The ion size is then gauged by the the fixed set of vdW parameters; the distance 
at which the ion-water vdW potential crosses zero minus the radius of a water molecule gives a rough indication of 
that size for the three anions. Coupled solute-solvent simulations were run for each of these solutes, and an additional 
coupled 500 ps sampling simulation including a Weeks-Chandler-Andersen (WCA) model potential^ 3 - with radius 
around 3 A was run for systems corresponding to the actual anions (1.1, 2.2, and 3.3). For comparison to experiment 
a sodium ion simulation (labeled by +) was also carried out. Uncharged, nonpolarizablc solutes corresponding to the 
above particle sizes were also simulated and have been denoted by the prefix m. See Table U for a compilation of the 
models studied in this paper. Properties of pure water were determined using reference simulations including WCA 
particles with radii of (pure water), 2.7, 3.1, and 3.4 A. In total, 20 systems were simulated requiring 21.2 ns of 
MD, after the exclusion of system 1.3 (Cl~-like ion with a = 7.25 A 3 ), which exhibited sporadic numerical divergence 
during the induced dipole calculations. 

The Amber implementation of AMOEBA was validated against the Tinker program and found to give identical 
energy components when Amber's long-range Lennard-Jones correction was turned off. Using Amber and systemati- 
cally varying the Ewald parameter, rj, between 0.4 and 0.6 A -1 for a chloride-water configuration carrying a net —1 
charge revealed a drift in the electrostatic energy on the order of 0.28 kcal/mol that could be removed by the addition 
of a uniform background charge correction term, 

E c = -332.0716-^- kcal/mol, (8) 

where V is the system volume and Q is its net charge. All energies reported here include this extra correction term 
(on the order of —0.4 to —0.5 kcal/mol depending on the system volume) not originally present in the raw Amber 
output. A discussion of this net charge correction is given in Rei£l (Eq. 5.50). 

For each system, the electrostatic potential at the center of the solute (labeled i) due to the surrounding waters 
was estimated using numerical differencing of E es (the Amber electrostatic energy plus E c ) for charges of +/ — 10~ 2 
electrons at zero solute polarizability. This potential was corrected for finite-size effects (Refs^l andS, Eq. 5.53) due 
to simulating with the Ewald potential at a solute charge q s i m via addition of g s i m £/L, where £ = —332.0716x2.837297 
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TABLE I: Summary of the systems modeled in the present study. The labels are defined in the text. 
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Finally, the anisotropy of the solvation environment around a given solute has been analyzed graphically by averaging 
the distance to the center of mass (COM) of the nearest n waters, generalizing a procedure used previously by other 
authors i 19 ' 68 First, all waters were rotated into a local coordinate frame defined using the COMs of the nearest 3 
waters. This was done in the usual way, using the first water COM to define the x axis and constraining the second to 
lie in the xy plane with positive y value. We reflect about the z axis to give the third water's COM a positive z value. 
Next, the total COM of the nearest n waters was averaged in this coordinate system. Since it was found that the 
major direction of change for this average with respect to n was in the < 1,1,1 > direction, all COMs were projected 
along this direction to create final one-dimensional anisotropy plots. These plots show the tendency of water to stack 
cither in the direction of the first three solvent waters (positive values) or on the opposite side (negative values). We 
have found that the anisotropy plots illustrate the variation in solvation structure more clearly than the unsigned 
distance distributions. 



IV. RESULTS 



A. Cavity formation free energies 



The AMOEBA water model has been parametrized to reproduce molecular multipole moments, heats of vapor- 
ization, and structures of small water clusters as well as bulk liquid density and internal energy i 40 i 41 It reproduces 
these quantities rather well, obtaining in addition good agreement with liquid phase radial distribution functions, heat 
capacity (21 +/- 5 cal/mol/K), average dipole moment in the liquid (2.78 D), and dielectric constant (82). For un- 
derstanding ion hydration in the context of quasi-chemical theory, we have computed the cavity formation probability 
profile (po) as a function of A, which yields the packing contribution to the free energy. 

Figure [2] shows the hard sphere chemical potential profile derived from the cavity occupancy probability for the 
AMOEBA and SPC/E water models. The major difference between the SPC/E and AMOEBA water models for 
these profiles is an upward shift in the AMOEBA results that can be attributed to an increase in the liquid phase 
density at 1 bar from 0.9556 in SPC/E to 1.0171 g/mL in AMOEBA. This AMOEBA value for the density is slightly 
larger than the previous report of 1.0004 g/mL,— a difference which may be due to the slightly longer time step, 
larger dipole convergence cutoff, or choice of temperature coupling algorithm for our calculations. Despite this, the 
difference in /igg at A = 3.5 A is only 0.5 kcal/mol. From these results, it would appear that fixed-charge water 
models can provide accurate estimates of nonpolar solute hydration free energies. 
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FIG. 2: Free energy cost for creating an oxygen unoccupied cavity of size A, /i§s HS (A), in NPT simulations of AMOEBA or 
SPC/E water models. 

B. Partial Molar Quantities 

We have calculated the hydration free energies for all solutes studied in this report and find excellent overall 
agreement with both previously reported values and experimental whole salt free energies, as reported in Table [Til 
The hydration free energies for the NaCl and NaBr pairs computed with the AMOEBA force field are within 0.2 
kcal/mol of the recently tabulated experimental results in Refs^ andZ^. The AMOEBA result for the Nal pair is 
within 0.5-1.0 kcal/mol of the experimental values. The error estimates were taken simply as the standard deviation 
over the several A values used to compute the free energies. As shown below in Figure [3 the errors come largely from 
the estimates of the (small magnitude) inner-shell and packing contributions, since the long-ranged part converges 
rapidly as a mean-field estimate. Thus it is likely the error estimates given are conservative and over-estimate the 
actual error in the total free energy calculations. 

The partial molar volumes of the AMOEBA ions computed by differencing average simulation volumes agree 
moderately well with experiment within their stated uncertainties. Despite the whole salt agreement, the single-ion 
partial molar volume data stand in contrast to Ref^i, giving Na + a positive value at 5.0 + / — 2.3 cm 3 mol -1 , and Cl~ 
a correspondingly smaller value. In the next section we show that the thermodynamic size decreases with increased 
polarizability, and note here that for system 1.2 (CI - size with a = 5.65 A 3 ), the computed partial molar volume 
becomes slightly negative. 

Because it is experimentally impossible to measure single-ion partial molar quantities, proton hydration is usually 
employed as a reference process. The present simulations can be used to infer values for this process. The infinite 
dilution free energies and enthalpies reported here, however, neglect issues associated with the average potential of 
the liquid phase because they lack a liquid/vapor interface. Energies calculated using this surfaceless approach have 
been termed 'intrinsic' and can be corrected for the presence of a surface potential by adding q<p w \. 

The proton solvation quantities were predicted by combining our data with experimental H + -Na + differences as 
well as direct measurements of halide acids. Intrinsic free energies calculated in this way vary little over estimates 
from alternative ions, giving —252 + </> v i kcal/mol with a simulation error of less than 1 kcal/mol, in relatively good 



8 





Na-Cl 


Na-Br 


Na-I 


Current 


-174.2 (1.6) 


-167.7 (2.1) 


-159.7 (2.0) 


Friedman^ 


-174.1 


-170.8 


-159.7 


Tissandies^ 


-174.0 


-167.6 


-158.7 


Schmidt 


-174.0 


-167.6 


-159.2 


dV/dN 


Current 


20.5 (3.3) 


22.5 (3.4) 


41.1 (3.2) 


Friedman 


16.6 


23.5 


35.0 


Krumgala^ 


16.62 


23.479 


34.998 



TABLE II: Partial molar hydration quantities for whole salts at infinite dilution, in kcal/mol and cm 3 mol 1 . Numbers in 
parentheses are one standard deviation errors. 



agreement with earlier work i 39 i 7 i 58 i 73 

The surface potential also creates a possible ambiguity in the partial molar volume. V = a g r ° al — ^ + q^p^, 
emphasizing the dependence of the surface properties on the experimental conditions. The data on partial molar 
volumes for HC1 and sodium salts have been used to calculate an intrinsic proton partial molar volume. The computed 
estimates made using these different experimental values, however, vary between 0.1 and 6.2 cm 3 mol -1 with our 
simulation uncertainty on the order of 2.3 cm 3 mol" 1 . 

C. Free Energy Components 

As discussed above, quasi-chemical partitioning has been used to separate the total hydration free energy into 
physical components. The ion-water radial distribution functions (rdfs) are shown in Figure [3] for reference. The 
A-depcndcnccs of the chemical contribution from inner-shell waters for each anion studied are shown in Figure 5] 
At first, it may be surprising that there is no identifiable jump due to breaking of direct-contact interactions with 
water which would be expected around 3.1-3.3, 3.3-3.45, or 3.5 — 3.7 A for CI", Br", or I - , respectively! 7 -^ The above 
distances correspond to locations just beyond the first minimum in the ion- water (hydrogen) rdf, or roughly at the first 
maximum of the ion- water (oxygen) rdf. Instead, the data show a smoothly increasing energetic penalty for pushing 
waters away from the anion - indicating the role of more distant waters (perhaps through many-body polarization 
effects) in creating a positive potential at the ion's center that compensates for the disruption of direct contacts. 

Increasing the ion polarizability results in a systematic downward shift in the /^s (A) profile, reflecting an increased 
probability of inner-shell water occupancy. This provides direct evidence for the decrease in thermodynamic ion size 
with increasing polarizability previously inferred by several authors ^ 19 ' 22 ! 50 This finding is also reflected by the trend 
in partial molar volumes discussed above. It is well known that the introduction of charge on a neutral solute has an 
even larger impact on the local solvation structure. 

The His(X) profiles presented in Figure|4]are computed out to A values that are just beyond the location of the first 
maximum in the ion- water (oxygen) rdfs. For the total free energy calculations discussed below, the A m f values used 
were taken to be the average distances to the closest water oxygen, (i? m in)i, which are roug hly 3.12, 3.29, and 3.5 A 
for the CI", Br", and I" ions, respectively. The magnitudes of /ijg for these values of A are quite small (on the order 
of -0.5 kcal/mol or less), and exhibit little polarizability dependence. 

Because the inner-shell contribution is so small for the chosen A m f (and shows little ion polarizability dependence) , 
and the packing term is a well-characterized size-dependent contribution independent of ion polarizability, the origin 
of the observed free energy differences between ions can be investigated by examination of the outer shell, long-ranged 
free energy component (Figure[5]). As described above, for A values large enough to yield mean-field behavior, ^ s LR 
is the average of the mean interaction energy between the coupled and uncoupled simulations. This allows for a simple, 
physically motivated division of the free energy for this step in terms of the first-order electrostatic (es), induction 
(ind), and van der Waals (vdW) energies. 

The division of AU is the usual one for force field simulations, with the exception that the electrostatic energy 
difference, AE CS = E cs (solvent + solute) — E cs (solvent), includes many-body induction due to the instantaneous 
changes in the system dipoles in the presence of the solute charge, qt. We define AE CS as 

AE CS = q&i + AE ind . (10) 

The potential at the solute charge location due to the surrounding solvent's self-consistent dipolc configuration (but 
uncoupled from the solute), <f>i, is as defined in Eq. [9j The Ai?i„d term includes the remainder of the electrostatic 
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FIG. 3: Ion-water(hydrogen) and ion-water(oxygen) radial distribution functions for the Na + , CI , Br , and I ions modeled 
with the AMOEBA force field. 



energy change, which is the energy change to repolarize the ion and the waters as they mutually interact. We note 
that a size-dependent correction (£qf/2L, where £ = —2.837297 and L is the box size) is omitted from AE es here 
since it has no impact on the computed free energies. 

For nonpolarizable force fields, the induction energy is zero. For polarizablc force fields, the potential/induction 
energy separation can be likened to the symmetry-adapted perturbation theory (SAPT) 75 i 76 ' 77 partitioning of in- 
teraction energies employed in quantum calculations. Along these lines, we term the q^i piece the 'first-order 
electrostatic energy' and the remainder will be called the induction contribution. The remaining contribution 
A-Eydw = -E v dw( s °l ven t + solute) — -E v dw (solvent) has the traditional interpretation of repulsion plus pairwise disper- 
sion energy. These three parts, then, make up the total solute-solvent interaction energy for any given configuration 
of the classical model: 



AU = qi<Pi + AE ind + AE- 



vdW- 



(11) 



Plots of the estimated free-energy contributions obtained from this partitioning as functions of the solute size and 
polarizability, a, are presented in Figure [3 The energy components have been calculated at a A corresponding to the 
average distance to the nearest water oxygen (J? m i n ) 1 as listed above. The free energies calculated from these radii 
are in good agreement with the estimates over all A used in creating Table HT1 

As can be seen in the upper-left panel of Figure the total free energy and total long-ranged contribution show 
a much larger dependence on ion size than on polarizability, consistent with free energy calculations using other 
polarizable models 4£ The addition of the /z|g + /Ltgg HS combination does not alter the polarizability dependence but 
does exhibit some specificity through a purely size-dependent effect in /iQg HS (above) , further separating the hydration 
free energy differences. The first-order electrostatic contribution is the largest part of the long-ranged contribution, and 
it exhibits size-dependent specificity and a slight decrease in favorability with increasing polarizability. As expected, 
the induction contribution displays the largest dependence on polarizability, and the observed downward trend with 
increasing polarizability overcomes the reverse trend in the first-order electrostatic part. The vdW piece exhibits an 
ordering that is inverted relative to the electrostatic parts, with the Cl~ ion displaying the largest positive value. 
Thus this vdW part of the long-ranged contribution serves to counteract somewhat the size-driven separation in free 
energies due to the electrostatic terms. It appears that the attractive dispersion component for the larger Br~ and 
I - ions counteracts the more repulsive contribution for the Cl~ ion. Interestingly, the combination of the outer-shell 
packing and the vdW part of the outer-shell, long-ranged contributions (/Xgg hs + ^os lr( v( 1W)) shows little specificity 



10 




2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 

X/A 

FIG. 4: Effect of polarizability on the energy of collapsing the inner-shell partition. For ions of all sizes (first index: 
l=Cl,2=Br,3=I), the curve shifts downward with increasing polarizability (second index: 0=zero,l=Cl,2=Br,3=I) indicating 
that increased polarization causes water to pack more closely. 

for the Cl~,Br~, I - sequence: 11.2, 11.1, and 11.9 kcal/mol, respectively. Since the also displays little specificity 
for the chosen A values, these results suggest that most of the specificity is carried in the electrostatic contributions 
to /LtQg LR in the classical model. 

In order to better understand the relative roles of the various parts of the electrostatic interactions, we here 
examine the electrostatic potentials at cavity and ion centers under various circumstances. The orientation of 
water molecules in the vicinity of a neutral cavity gives rise to a positive electrostatic potential at its center i 67 ' 78 This 
potential is plotted as a function of A, the distance to the nearest water oxygen, in the lower section of Figure [5] 
The AMOEBA potential is slightly higher (« 60 mV) than that at the center of an uncharged solute^ 8 - because our 
cavity definition allows intrusion of water hydrogens, which would otherwise be repelled. Although both AMOEBA 
and SPC/E water models share a similar trend, the charge distributions in the AMOEBA model serve to shift the 
cavity potential downwards by sa 100 mV. This is interesting because, while the AMOEBA fixed partial charges are 
reduced in magnitude relative to the SPC/E model, the average dipole magnitudes for the AMOEBA waters near the 
cavity (2.76 D) are close to the AMOEBA bulk value (2.78 D) and are larger in magnitude than the fixed permanent 
dipole of the SPC/E model (2.35 D). These AMOEBA dipole magnitudes near the cavity are in sharp contrast to the 
marked decrease in the dipole moments approaching the water/vacuum interfaced We can conclude that many-body 
polarization effects are important in orienting the dipolcs around small cavities. 

The first-order electrostatic contribution to the long-ranged part of the free energy comes from the average of 
the potential for the uncoupled and coupled cases. To examine the coupled case, we have added to Figure [5] the 
electrostatic potential at the center of ionic solutes (Na + and Cl~ ) as a function of the distance to the nearest water 
oxygen. The ion models used with the SPC/E water model were the same as in our previous publication^ while the 
AMOEBA ions are taken directly from the present study. First, the first-order potential <I>i for the coupled case is 
much larger for the SPC/E model than for the AMOEBA model, reaching values close to -9 V for the Na + ion at 
small A values. Second, there is a clear break in slope for the SPC/E model at A values close to the first maximum 
in the ion- water (oxygen) rdf; pushing the waters further away from the favorable solvation environment creates a 
large energetic penalty. The polarizable AMOEBA model, on the other hand, softens this effect through many-body 
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FIG. 5: Contributions to the OS-LR conditional free energy average at A ~ (-Rmin)i- As explained in the text, the components 
are averages over the coupled and uncoupled cases. The corresponding contributions from a nonpolarizable C1~-SPC/E system 
are shown as single points. In the upper left-hand figure, the upper curves for each ion are the total hydration free energies 
obtained by adding the inner-shell and outer-shell packing contributions. 



polarization effects. For the SPC/E case, the total electrostatic part of the long-ranged free energy contribution 
comes entirely from the first-order electrostatic contribution. The AMOEBA model makes up the remainder of the 
free energy through induction effects, and these are seen to be substantial in magnitude (about 25% of the total 
electrostatic part of the free energy). Making the above plot with a nonpolarizable CI - ion in AMOEBA does not 
significantly change the shape of the potential profile, but shifts it upward by « 140 mV. We conclude that many- 
body polarization contributions are significant in representing the electrostatic potential in nonpolar cavities and ion 
centers. The cavity free energy is rather insensitive to polarization, however, as seen in Figure [5] 

As a point of interest, the complete electrostatic contribution to the long-ranged part of the free energy can be 

2 

compared to the Born model, jugg LR os ps qi&o(\) — (1 — 1/e) J# (corrected for nonzero cavity potential, ^oW, Figure 
[5]); R is A — Rwj where R\y is the radius of the water molecule, taken here as 1.3 A. When such a comparison is 
made (Figure [7]), we see that in a general sense the Born model is relatively accurate even at these short length scales 
near the ion. Since the Born model is equivalent to the Gaussian model employed here for the OS,LR contribution^ 
the agreement is not surprising. On the other hand, it is clear that molecular details matter in obtaining an accurate 
representation of the electrostatic part of the free energ y 79 ' 80 ! 81 since the slopes of the AMOEBA-generated curves 
differ from the Born model, and there are noticeable jumps between ions for a given A value. 



D. Analysis of local solvation structure 



In addition to the thermodynamic and electrostatic analysis above, it is interesting to explore the local solvation 
structure around the various ions as functions of size and polarizability. In Figure [8l the solvation-shell water ordering 
is addressed by plotting the signed center of mass distance of solvating water clusters of increasing size (above), 
ordering waters using the distance from the ion to their center of mass. A recent report^ examined the same question 
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by projecting along the direction of the instantaneous ion dipole. Yet another possibility would be to search directly 
for molecular cavities in the solvation water by, for example, computing a potential of mean force between the ion 
and a hard sphere. 

This anisotropy is significant because it shows one of the mechanisms by which ordering may increase in response 
to increased polarization. It also indicates a possible role of polarization on surface affinity, since a rightward shift of 
the first zero in the anisotropy plot may indicate a higher solute preference for surface-like solvation shells^ 

We measure the degree of solvation asymmetry by the number of waters just before the anisotropy plot reaches 
its first zero. Table IIIII presents the numerical values for these cluster sizes, along with a range of other solvation 
data. Figure [5] displays the solvation anisotropy for a range of cases. In the upper-left panel, it is apparent that the 
water molecule and the Na + ion exhibit relatively symmetric solvation shells in comparison with the three anions and 
nonpolar solutes. For the three anions, it is clear that increasing the ion polarizability shifts the curves rightward, 
indicating increased anisotropy. The anisotropy also increases through the sequence CI - , Br~, I™, for their normal 
AMOEBA parameter values (systems 1.1, 2.2, and 3.3). In the lower-left panel, we present the anisotropy plot for the 
sequence CI™, Br™, I™ but with a fixed ion polarizability of a = 5.65 A 3 . Surprisingly, there appears to be little size 
dependence to the anisotropy. The final lower-right panel displays the anisotropy for uncharged solutes with vdW 
parameters corresponding to each ion. For uncharged solutes, there is a clear size dependence, and the water solvation 
shell is quite asymmetric. 

Table Hill collects the solvation environment properties relevant to this discussion. The total hydration free energies 
are also presented for each case for reference. The number of waters just before the anisotropy plot first crosses zero, 
Ao, is largest for the uncharged, non-polarizable solutes indicating the highest asymmetry for those cases. The Na + 
ion exhibits the least asymmetry, and No appears to correlate most strongly with polarizability for the anion series. 
The most probable solvation number (the number of waters with hydrogens within 3.8 A of the ion and with an 
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sizes produce about a 5 kcal/mol separation when compared at equal nominal radii, contradicting the Born model prediction. 
Labels are as in Figure [4] 
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TABLE III: Comparison of local solvation environment indicators: No, the number of waters required for Figure [8] to reach 
zero, average distance to the closest water oxygen (A), average ion dipole moment (Debye), average first-shell water dipole 
moment (Debye), and single-ion total hydration free energy (kcal/mol). Numbers in parentheses indicate numerical uncertainty 
in the last digit. The label m indicates an uncharged solute. 
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ion-H-0 angle greater than 130°) is 6 for all of the anion cases (data not in table). The average distance to the closest 
water oxygen shows relatively little dependence on ion polarizability, but does show a substantial compression from 
turning on the ion charge. The average ion dipole moment magnitude exhibits strong dependence on ion polarizability 
as expected. As points of reference, the average dipoles for the Cl _ and Br~ ions with their AMOEBA parameters are 
1.56 and 2.19 D, respectively, compared with the AIMD values of 0.8 and 1.0 D^i It thus appears that the AMOEBA 
force field is appreciably over-polarizing the ions. The average dipole magnitudes for the waters in the first solvation 
shell (criterion given above) are slightly suppressed compared to bulk water (2.79 D). Recent AIMD studies have 
made similar observations for both anions^ and cations^ 



V. SUMMARY 

This paper has presented an analysis of the thermodynamics and local solvation structure involved in polarizable 
anion hydration. The theoretical tool employed for the analysis is the quasi-chemical theory of molecular solutions.— 
The present analysis shows that, in these classical polarizable models, the predominant factors affecting hydration 
free energies are the ion charge, size, and polarizability, in decreasing order of importance. Varying the ion size and 
polarizability independently is physically artificial, but serves to elucidate the relative importance of these factors. 
The QCT partitions the free energy into three components: inner shell, outer-shell packing, and outer-shell long- 
ranged contributions.— That division is exploited here by noting that, as the conditioning radius is increased, the 
outer-shell long-ranged contribution assumes a Gaussian character ! 60 ! 63 The Gaussian nature of the coupling energy 
distributions in turn implies that a mean-field treatment should be accurate. This leads to a natural division of the 
long-ranged component of the hydration free energy into first-order electrostatic, induction, and van der Waals parts. 

We have found that, for the sodium and halide ions examined here, the conditioning radii can be taken as quite 
small. This is likely due to the spherical nature of these simple ions; future work will focus on application of this 
approach to molecular ions such as the nitrate ion. With these small conditioning radii, the /iff + /igg HS (inner shell 
plus packing) combination is small and dominated by the packing part. Thus fif^ + /i™ s HS exhibits ion specificity 
mainly through the purely size-dependent packing piece. We also found that the sum of the inner-shell, packing, and 
long-ranged vdW parts of the free energy is roughly constant, suggesting electrostatic effects (induced by ion size) 
are the primary determinants of specificity for the chosen hard-sphere conditioning radii. The inner-shell piece does 
increase somewhat in magnitude (becomes more negative) with increasing ion polarizability, indicating electrostriction 
of nearby waters, but this effect is small for the conditioning radii used here. The main contributions to the hydration 
free energy come from the first-order electrostatic and induction terms. The induction contribution comprises about 
25% of the electrostatic part of the free energy. The first-order electrostatic piece becomes slightly less favorable with 
increasing ion polarizability, while the induction piece becomes more favorable. Due to the induction contribution, 
the total free energy becomes slightly more favorable with increasing ion polarizability. 

Comparison of the electrostatic potential at the center of nonpolar cavities and ions illustrates substantial differences 
between the fixed-charge SPC/E and polarizable AMOEBA models. For the nonpolar cavity case, the electrostatic 
potential at the cavity center is reduced in the polarizable AMOEBA water model compared to the SPC/E water 
model. This occurs even though the average total dipole magnitudes for the nearby waters are larger for the AMOEBA 
model than for the SPC/E model. This illustrates the importance of many-body effects in forming the electrostatic 
environment near the cavity. Those differences are accentuated when the average potentials arc examined for the Na + 
and Cl~ ions. In addition, the SPC/E model exhibits a clear break in the slope of the potential vs. A plot at radii 
close to the first maximum of the ion- water (oxygen) rdf, while the polarizable model docs not display such a break. 
Apparently the many-body polarization smooths out the energy changes as waters are pushed away from the ions. To 
examine further the electrostatic contribution to the hydration free energy, we compared the long-ranged electrostatic 
contribution to the Born model of solvation. Those results show that, while the Born model fits the simulation data 
well considering its simplicity, it fails to capture molecular-level details important for accurate computations of free 
energies. 

Finally, we found that, contrary to the hydration free energies, the local solvation structure is affected more by 
ion polarizability than by size. Increasing the ion polarizability leads to more asymmetric solvation environments. 
Changes in size for a given polarizability lead to quite small changes in the solvation anisotropy. The average dipole 
moment magnitudes for waters in the first solvation shell are slightly suppressed from the bulk value. This highlights 
the importance of many-body water-water interactions in responding to the strong field of the nearby ion. The average 
ion dipoles, however, displayed values considerably larger (by a factor of 2) than values observed in AIMD quantum 
simulations^ 
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VI. DISCUSSION 

The results presented here for the hydration free energies show good agreement with experimen t 69 ' 70 and previous 
simulations employing the AMOEBA models Based on the free energy calculations alone, the AMOEBA model 
appears to provide an excellent representation of the hydration of simple alkali halidc ion pairs. The quasi-chemical 
approach allows us to go further, however, in disentangling the various contributions to the free energies. In partic- 
ular, the quasi-chemical partitioning has allowed us to examine in detail the relative contributions of ion size and 
polarizability in hydration. 

Examination of the structural properties of the surrounding solvation shell points out some potential limitations 
of the AMOEBA model, and likely other models which possess limited or no damping of close-ranged electrostatic 
interactions. Namely, as discussed above, the anions are over-polarized compared with AIMD simulations , 54 ' 68 ' 83 other 
quantum models,— and recent carefully parametrized classical models.— In another paper supporting this point, 
we have examined the local structure in related QM/MM calculations^ 4 - In that work we generated an ensemble 
of configurations from an AMOEBA simulation and then performed MP2-levcl electronic structure calculations on 
clusters extracted from the ensemble. The ion plus the nearest 6 waters were treated quantum mechanically, while 
the more distant waters produced an electrostatic potential due to their AMOEBA-generated multipole moments. 
The quantum charge distributions were analyzed with both the ChelpG and QTAIM methods^ We found that, in 
agreement with the AMOEBA simulations, the closest waters possess dipoles close to the bulk value due to collective 
interactions between the water molecules. The anion dipoles are much smaller (roughly 0.6 D for the Cl _ ion) than 
those predicted by the AMOEBA force field (1.56 D), however. In addition, we observed substantial charge transfer 
(of magnitude 0.2 charges) from the anion to the surrounding waters. This finding is consistent with several previous 
papers , 86 ' 87 ' 88 ' 89 but differs from the conclusions of Ref.— that applied a different charge-accounting scheme. We 
note that the anion dipole estimate in Reb^ 4 - is somewhat smaller than the AIMD prediction; this may be due to 
the observed charge transfer which is not accounted for in the Wannier analysis typically employed with AIMD 
simulations^ 4 - 

The majority of classical simulations that have examined surface affinity of anions have targeted the ion polarizability 
as the key determinant of the observed segregation^^ 1 - Other simulations have emphasized the important role of ion 
size^i 5 - while Hagberg, Brdarski, and Karlstrom 7 - have shown that both factors play important roles. Herce et ai.jSP. 
using polarizable models, suggested that the main causes of surface solvation of a single ion in a water cluster are 
both water and ion polarization, coupled to the charge and size of the ion. On the other hand, Wick and Xantheas 5 - 5 - 
have linked anisotropic local solvation structure to surface affinity - perhaps the ions exhibiting the most asymmetric 
solvation shells prefer the highly asymmetric liquid- vapor interface region. Intuitively, many-body polarization should 
be a major contributor ; 55 ' 90 especially when the large electric field in the region of the water/vapor interface is noted. 
The observed over-polarization in some classical models and charge transfer in quantum studies both suggest that 
further improvements in the representation of polarization in classical models should be considered to help to quantify 
the forces driving surface affinity of anions. A recent study by Wick 9 -^ has confirmed that increased damping of the 
polarization interactions, leading to reduced ion polarization (in agreement with AIMD simulations), in fact results 
in reduced anion surface affinity. Wick's work is consistent with previous results by Ishiyama and Morita— that 
employed reduced polarizability models. 

A first remedy for the anion over-polarization in the classical models would be to reduce the anion polarizability 
employed until agreement is obtained with ab initio results for the average dipole magnitudes. Masia 5 - 2 - has argued 
against this approach, however, since this does not properly represent the physics of the dipoles at close contact; a 
large magnitude tail can still exist in the distribution of dipole moments due to over-polarization from closeby waters. 
A second remedy would be to re-parametrize the Thole-type damping functions to mimic more closely the longer- 
ranged damping functions developed by Masia.— A third remedy would be to focus on a more accurate QM/MM 
approach as outlined above^ 4 - We are currently developing such a QM/MM (MP2 level) approach for computing 
the free energy contributions in the context of quasi-chemical theory. We have chosen this level of theory so that 
many-body dispersion interactions as well as accurate representations of local charge rearrangements can be included 
accurately. Such studies should help to elucidate further the driving forces for ion surface segregation, and allow for 
tests of recently developed theorie a 24 ' 26 of this phenomenon. 
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